// --------------------------------------------------------------------
//
//  File:        OAk8D2.rtf [Optimal Ate pairing for twisted Edwards curves with quadratic twists]
//  Date:        Oct. 10, 2018
//  Last update: Oct. 10, 2018
//  Description: Magma implementation of the optimal ate pairing on twisted Edwards curves with quadratic twists.
//               Embedding degree k = 12, CM discriminant D = 2
//
//  (c) 2018, Georgios Fotiadis & Chloe Martindale
//                 
//  Georgios Fotiadis, University of the Aegean, Greece, email: gfotiadis@aegean.gr
//  Chloe Martindale, Technische Universiteit Eindhoven, the Netherlands, email: chloemartindale@gmail.com
// --------------------------------------------------------------------
//  Polynomial Family
//  p(x) = (x^14 − 4*x^10 + 2*x^8 + 4*x^6 − 2*x^4 + 5*x^2 + 2)/8
//  t(x) = x^2 + 1
//  r(x) = x^8 - x^4 + 1 
// --------------------------------------------------------------------
//  Optimal ate pairing: a(Q,P) = [ f_{x^2,Q}(P) ]^((p^k - 1)/q)
// --------------------------------------------------------------------

//Weierstrass to Montgomery
MontPoint:=function(P, alpha, s)
  x:=s*(P[1] - alpha);
  y:=s*P[2];
  MP:=Vector(2,[x,y]);
  return MP;
end function;

//Montgomery to Edwards
EdPoint:=function(P)
  X:=(P[1]/P[2]);
  Y:=(P[1] - 1)/(P[1] + 1);
  Z:= 1 ;
  T:=(X*Y)/Z;
  EP:=Vector(4,[X,Y,Z,T]);
  return EP;
end function;

//Miller's function point doubling: h_{P,P}(Q) 
hDBL:=function(P, Eda, Edd, eta, y0)
  X1:=P[1]; Y1:=P[2]; Z1:=P[3]; T1:=P[4];
  A:=X1^2;
  B:=Y1^2;
  C:=Z1^2;
  D:=(X1 + Y1)^2;
  E:=(Y1 + Z1)^2;
  F:=D - (A + B);
  G:=E - (B + C);
  H:=Eda*A;
  I:=H + B;
  J:=C - I;
  K:=J + C;
  X3:=F*K;
  Y3:=I*(B - H);
  Z3:=I*K;
  T3:=F*(B - H);
  CZ2:=2*Y1*(T1 - X1);
  CXY:=2*J + G;
  CXZ:=2*(Eda*X1*T1 - B);
  R:=Vector(4,[X3,Y3,Z3,T3]); 
  h:=CZ2*eta + CXY*y0 + CXZ;
  return R, h;
end function;

//Miller's function point addition h_{P,S}(Q) 
hADD:=function(P, S, Eda, Edd, eta, y0)
  X1:=P[1]; Y1:=P[2]; Z1:=P[3]; T1:=P[4];
  X2:=S[1]; Y2:=S[2]; Z2:=S[3]; T2:=S[4];
  A:=X1*X2;
  B:=Y1*Y2;
  C:=Z1*T2;
  D:=Z2*T1;
  E:=D + C;
  F:=(X1 - Y1)*(X2 + Y2) + B - A;
  G:=B + Eda*A;
  H:=D - C;
  I:=T1*T2;
  X3:=E*F;
  Y3:=G*H;
  Z3:=F*G;
  T3:=E*H;
  CZ2:=(T1 - X1)*(T2 + X2) - I + A;
  CXY:=X1*Z2 - X2*Z1 + F;
  CXZ:=(Y1 - T1)*(Y2 + T2) - B + I - H;
  R:=Vector(4,[X3,Y3,Z3,T3]); 
  h:=CZ2*eta + CXY*y0 + CXZ;
  return R, h;
end function;

//final exponentiation k = 12, D = 2
finalexp:=function(f, x0, p, q, Fp)
  t0:=1;
  f:=Frobenius(f, Fp, 6)/f; 
  f:=Frobenius(f, Fp, 2)*f; 
  x2:=x0^2;
  x3:=x2^2;
  x6:=x2^3;
  l:=(x6 - 3*x2 + 2) div 8;
  y:=f^(l);
  t:=Frobenius(y, Fp, 3);
  t0:=t0*t;
  t:=Frobenius(y^x2, Fp, 2);
  t0:=t0*t;
  t:=Frobenius(y^(x3 - 1), Fp, 1);
  t0:=t0*t*y^(x6 - x2)*f;
  return t0;
end function; 

//Optimal ate pairing calculation
OptAtePairing:=function(Q, V, q, d, m, p, k, z0, eta, y0, Fp)
  f:=1;    
  T:=Q;
  s:=Intseq(Abs(V[1]),2);
  for i:=#s - 2 to 0 by -1 do
    T,h:=hDBL(T, d, m, eta, y0); 
    f:=h*f^2; 
    if s[i+1] eq 1 then
      T,h:=hADD(T, Q, d, m, eta, y0); 
      f:=h*f; 
    end if;
  end for;
  f:=finalexp(f, z0, p, q, Fp);
  return f;
end function;

//MAIN PROGRAM
t2:= 0 ; LOOP:= 100 ;
D := 2 ; k := 12 ;                 
p := 90898520127544074556538387695336173057886717639583659376403912072593778341743710610181181036677380520056745403873138441122212347150337 ;
q := 115823896522048476472095240657222785056166991332278538593381683991139536404481 ;      
h := 784799362282207836808044004755169147455747243193650905088 ;    
z0:= 4295114753 ; 
a:= 10353841140864156393850331915848774030682220612200405154255017002541059163488608120808766441757343152957214120074465669038166665814946 ; 
b:= 39473289363775963104102109451594297948691152843912940945273433085571876583864616864154630093167743249959592989117035943541619831594138 ;
V:=Vector(4,[z0^2,-1,0,0]) ; //the shortest vector of the lattice L

s:=k div 2;
Fp := GF(p);
P<x> := PolynomialRing(Fp);

E:=EllipticCurve([Fp!(a), Fp!(b)]); 
f:=x^3 + a*x + b; 
n:=h*q; 
t:=p + 1 - n; 

ROOTS:=Roots(f);
ALPHA:=Vector(3,[ ROOTS[1][1],ROOTS[2][1],ROOTS[3][1] ]);

for i:= 1 to 3 do
  test,sqrt:=IsSquare(1/(3*ALPHA[i]^2 + a));
  if (test eq true) then
    ss:=sqrt;
    alpha:=ALPHA[i];  
    break i;
  end if;
end for;

//Montgomery Curve
A:=3*alpha*ss; 
B:=ss;
//Twisted Edwards Curve
Eda:=Fp!((A + 2)/B); 
Edd:=Fp!((A - 2)/B);

P<z> := PolynomialRing(Fp);
g:=IrreduciblePolynomial(Fp,s); 
Fps<u> := ExtensionField<Fp,z|g>;
R<zz> := PolynomialRing(Fps);

for i1:= 1 to 5 do
  for i0:= 0 to 5 do
    fX:=i1*u + i0; 
    FX:=zz^(k div s) - fX;
    if IsIrreducible(FX) eq true then
      Fpk<v> := ExtensionField<Fps,zz|FX>; 
      x0:=Roots(FX, Fpk); 
      w:=x0[1][1];
      break i1;
    end if;
  end for;
end for;

Et:=EllipticCurve([Fps!(a*w^4),Fps!(b*w^6)]);
ht:=#Et div q;

alphat:=Fps!(alpha*w^2); 
sst:=Fps!(ss/w^2);
//Montgomery Curve
At:=3*alphat*sst; 
Bt:=sst;
//Twisted Edwards Curve
Edat:=Fps!((At + 2)/Bt); 
Eddt:=Fps!((At - 2)/Bt);

Ek:=EllipticCurve([Fpk!(a), Fpk!(b)]); 
Zk<xx> := PolynomialRing(Integers());
Nk := Resultant(xx^k - 1, xx^2 - t*xx + p);
hk := Nk div q^2;

Edak:=Fpk!(Eda); 
Eddk:=Fpk!(Edd); 
alphak:=Fpk!(alpha); 
ssk:=Fpk!(ss); 

for i:= 1 to LOOP do
//create two points PW and Pa in E(Fp) of order q (two points because of the bilinearity check)
G:=Random(E); 
nP:=h; 
P:=nP*G; 
nP:=Random(1, q - 1); 
Pa:=P; 
PW:=nP*P; 

MP:=MontPoint(PW, alpha, ss); 
EP:=EdPoint(MP);
MPa:=MontPoint(Pa, alpha, ss); 
EPa:=EdPoint(MPa);

eta:=w*(EP[3] + EP[2])/(EP[1]);
y0:=(EP[2]/EP[3]);

etaa:=w*(EPa[3] + EPa[2])/(EPa[1]);
y0a:=(EPa[2]/EPa[3]);

//create two points QW and Qa in Et(Fp^6) of order q (two points because of the bilinearity check)
G:=Random(Et); 
nQ:= ht ; 
Q:=nQ*G;  
QW:=Q; 
xQ:=QW[1]*w^(-2); 
yQ:=QW[2]*w^(-3); 
E1:=Ek![xQ,yQ]; 
XQ:=Frobenius(xQ, Fp, 1); 
YQ:=Frobenius(yQ, Fp, 1); 
D1:=Ek![XQ,YQ]; 
QW:=D1-E1; 

xQ:=QW[1]*w^(2);
yQ:=QW[2]*w^(3);
QW:=Et![xQ,yQ]; 

nQ:=Random(1, q - 1); 
Qa:=QW; 
QW:=nQ*QW; 

MQ:=MontPoint(QW, alphat, sst); 
EQ:=EdPoint(MQ); 

MQa:=MontPoint(Qa, alphat, sst); 
EQa:=EdPoint(MQa); 

nQP:=nQ*nP;
t1 := Cputime();
opt:=OptAtePairing(EQ, V, q, Edat, Eddt, p, k, z0, eta, y0, Fp); 
t2:=t2 + Cputime(t1); 
end for;
t2:=t2/LOOP;
printf"average time for %o pairings = %o\n", LOOP, t2;

//Bilinearity check
OptAtePairing(EQ, V, q, Edat, Eddt, p, k, z0, eta, y0, Fp) eq OptAtePairing(EQa, V, q, Edat, Eddt, p, k, z0, etaa, y0a, Fp)^(nQP);